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JONATHAN J. CROFTS AND RUSLAN L. DAVIDCHACK t 

Abstract. An algorithm for detecting periodic orbits in chaotic systems [Phys. Rev. E, 60 
(1999), pp. 6172—6175], which combines the set of stabiUsing transformations proposed by Schmelcher 
and Diakonos [Phys. Rev. Lett., 78 (1997), pp. 4733—4736] with a modified semi-implicit Euler it- 
erative scheme and seeding with periodic orbits of neighbouring periods, has been shown to be 
highly efficient when applied to low-dimensional systems. The difficulty in applying the algorithm 
to higher-dimensional systems is mainly due to the fact that the number of the stabilising trans- 
formations grows extremely fast with increasing system dimension. Here we analyse the properties 
of stabilising transformations and propose an alternative approach for constructing a smaller set of 
transformations. The performance of the new approach is illustrated on the four-dimensional kicked 
double rotor map and the six-dimensional system of three coupled Henon maps. 
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1. Introduction. Periodic orbits play an important role in the analysis of vari- 
ous types of dynamical system. In systems with chaotic behaviour, unstable periodic 
orbits (UPOs) form a "skeleton" for chaotic trajectories 3 . A well regarded definition 
of chaos [3 requires the existence of an infinite number of UPOs that are dense in 
the chaotic set. Different geometric and dynamical properties of chaotic sets, such 
as natural measure, Lyapunov exponents, fractal dimensions, entropies tl9j , can be 
determined from the location and stability properties of the embedded UPOs. Pe- 
riodic orbits are central to the understanding of quantum-mechanical properties of 
nonseparable systems: the energy level density of such systems can be expressed in 
a semiclassical approximation as a sum over the UPOs of the corresponding classi- 
cal system 9 . Topological description of a chaotic attractor also benefits from the 
knowledge of periodic orbits. For example, a large set of periodic orbits is highly 
constraining to the symbolic dynamics and can be used to extract the location of a 
generating partition [S] I22| . The significance of periodic orbits for the experimental 
study of dynamical systems has been demonstrated in a wide variety of systems |16| , 
especially for the purpose of controlling chaotic dynamics i20i| with possible applica- 
tion in communication jjS]. 

It is therefore not surprising that much effort has been put into the development 
of methods for locating periodic solutions in different types of dynamical systems. 
In a limited number of cases, this can be achieved due to the special structure of 
the systems. Examples include the Biham-Wenzel method applicable to Henon-like 
maps U, or systems with known and well ordered symbolic dynamics 11 . For generic 
systems, however, most methods described in the literature use some type of an 
iterative scheme that, given an initial condition (seed), converges to a periodic orbit of 
the chaotic system. In order to locate all UPOs with a given period p, the convergence 
basin of each orbit for the chosen iterative scheme must contain at least one seed. 
The seeds are often chosen either at random from within the region of interest, from a 
regular grid, or from a chaotic trajectory with or without close recurrences. Typically, 
the iterative scheme is chosen from one of the "globally" convergent methods of quasi- 
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Newton or secant type. However, experience suggests that even the most sophisticated 
methods of this type suffer from a common problem: with increasing period, the basin 
size of the UPOs becomes so small that placing a seed within the basin with one of 
the above listed seeding schemes is practically impossible |18| . 

A different approach, which appears to effectively deal with the problem of re- 
duced basin sizes has been proposed by Schmelcher and Diakonos (SD) |25[ I2f)| . The 
basic idea is to transform the dynamical system in such a way that the UPOs of the 
original system become stable and can be located by simply following the evolution 
of the transformed dynamical system. That is, to locate period-p orbits of a discrete 
dynamical system 

(1.1) U: /: M" R" , 
one considers an associated flow 

(1.2) E: ^ = Cgix), 

where g{x) — f^{x) — x and C is an n x n constant orthogonal matrix. It is easy to 
see that map f^{x) and flow E have identical sets of fixed points for any C, while C 
can be chosen such that unstable period-p orbits of U become stable fixed points of 
E. 

Since it is not generally possible to choose a single matrix C that would stabilise 
all UPOs of U, the idea is to find the smallest possible set of matrices C = {Ck}^^i, 
such that, for each UPO of U, there is at least one matrix C G C that transforms 
the unstable orbit of U into a stable fixed point of E. To this end, Schmelcher and 
Diakonos have put forward the following conjecture p5| : 

Conjecture 1. Let Csd be the set of all n x n orthogonal matrices with only 
±1 non-zero entries. Then, for any n x n non-singular real matrix G, there exists a 
matrix C G Csd such that all eigenvalues of the product CG have negative real parts. 

Observation 1. The set Csd forms a group isomorphic to the Weyl group 
Bn 114^ , i.e. the symmetry group of an n- dimensional hypercube. The number of 
matrices in Csd is K — 2"n!. 

The above conjecture has been verified for n < 2 and appears to be true for 
n > 2, but, thus far, no proof has been presented. According to this conjecture, any 
periodic orbit, whose stability matrix does not have eigenvalues equal to one, can 
be transformed into a stable fixed point of E with G G Csd- In practice, to locate 
periodic orbits of the map C/, we try to integrate the fiow E from a given initial 
condition (seed) using different matrices from the set Csd- Some of the resulting 
trajectories will converge to fixed points, while others will fail to do so, either leaving 
the region of interest or failing to converge within a specified number of steps. 

The main advantage of the SD approach is that the convergence basins of the 
stabilised UPOs appear to be much larger than the basins produced by other iterative 
schemes |2tilll5l lH]. making it much easier to select a useful seed. Moreover, depending 
on the choice of the stabilising transformation, the SD method may converge to several 
different UPOs from the same seed. 

The flow E can be integrated by any off-the-shelf numerical integrator. Schmelcher 
and Diakonos have enjoyed considerable success using a simple Euler method. How- 
ever, the choice of the integrator for this problem is governed by considerations very 
different from those typically used to construct an ODE solver. Indeed, to locate 
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a fixed point of the flow, it may not be very efficient to follow the flow with some 
prescribed accuracy. Therefore, local error considerations, for example, are not as 
important. Instead, the goal is to have a solver that can reach the fixed point in as 
few integration steps as possible. In fact, as shown by Davidchack and Lai j^, the 
efficiency of the method can be improved dramatically when the solver is constructed 
specifically with the above goal in mind. In particular, recognizing the typical stiff- 
ness of the flow E, Davidchack and Lai have proposed a modified semi- implicit Euler 
method: 

(1.3) = + - G,]-ig(x,) , 

where /3 > is a scalar parameter, Sj = \ \g{xj)\\ is an L2 norm, Gj = Dg{xj) is the 
Jacobian matrix, and "T" denotes transpose. Note that, away from the root of the 
above iterative scheme is a semi-implicit Euler method with step size h = {Psj)~^ 
and, therefore, can follow the flow E with a much larger step size than an explicit 
integrator (e.g. Euler or Runge-Kutta). Close to the root, the proposed scheme can 
be shown to converge quadratically 15 , analogous to the Newton-Raphson method. 

Another important ingredient of the algorithm presented in 0] is the seeding 
with already detected periodic orbits of neighbouring periods. This seeding scheme 
appears to be superior to the typically employed schemes and enables fast detection of 
plausibly all^ periodic orbits of increasingly larger periods in generic low-dimensional 
chaotic systems. For example, for the Ikeda map at traditional parameter values, the 
algorithm presented in ^ was able to locate plausibly all periodic orbits up to period 
22 for a total of over 10^ orbit points. Obtaining a comparable result with generally 
employed techniques requires an estimated 10^ larger computational effort. 

While the stabilisation approach is straightforward for relatively low-dimensional 
systems, direct application to higher-dimensional systems is much less efficient due 
to the rapid growth of the number of matrices in Csb- Even though it appears that, 
in practice, far fewer transformations are required to stabilise all periodic orbits of 
a given chaotic system jSTJ, the sufficient subset of transformations is not known a 
priori. It is also clear that the route of constructing a universal set of transformations 
is unlikely to yield substantial reduction in the number of such transformations. For 
instance, a smaller set of universal transformations with K — (n + 1)!, which is 
isomorphic to the Weyl group An, is sufficient to stabilise all types of periodic orbits 
for n < 4, but can be shown to fail for certain types of orbits when n> A. Therefore, a 
more promising way of using stabilising transformations for locating periodic orbits in 
high-dimensional systems is to design such transformations based on the information 
about the properties of the system under investigation. 

In this Article, we propose to construct stabilising transformations based on the 
knowledge of the stability matrices of already detected periodic orbits (used as seeds) . 
The advantage of our approach is in a substantial reduction of the number of trans- 
formations, which increases the efficiency of the detection algorithm, especially in the 

-"^It is not possible to prove, within our approach, the completeness of the detected sets of UPOs. 
Rather, our assertion of completeness is based on the plausibility argument. The following three 
criteria are used for the validation of the argument: 

i) Methods based on rigorous numerics (e.g. in have located the same UPOs in cases where such 
comparison is possible (usually for low periods, since these methods are less efficient). 

ii) Our search strategy scales with the period p (see J3]and |^). If we can tune it to locate all 
UPOs for low periods (where we can verify the completeness using (i)), it is likely (but not provably) 
capable of locating all UPOs of higher periods as well. 

iii) For maps with symmetries, we test the completeness by verifying that all the symmetric partners 
for all detected UPOs have been found (see H4.1l and i|4.2i . 
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case of higher-dimensional systems. The layout of the paper is as follows. In f]2we 
study the properties of the stabilising transformations for n — 2 and their relation- 
ship to the properties of the stability matrix of a periodic orbit. In Sj3|we extend the 
analysis to higher-dimensional systems and show how to construct stabilising trans- 
formations using the knowledge of the stability matrices of already detected periodic 
orbit points. In particular, we argue that the stabilising transformations depend es- 
sentially on the signs of unstable eigenvalues and the directions of the corresponding 
eigenvectors of the stability matrices. Section ^illustrates the application of the new 
stabilising transformations to the detection of periodic orbits in a four-dimcnsional 
kicked double rotor map and a six-dimensional system of three coupled Hcnon maps. 
We conclude with the summary and discussion of possible further developments of 
the stabilising transformations approach in ^ 

2. Stabilising transformations in two dimensions. The stability of a fixed 
point X* of the flow S is determined by the real parts of the eigenvalues of the matrix 
CG, where G = Dg{x*) is the Jacobian matrix of g{x) evaluated at x* . For x* to be a 
stable fixed point of S, the matrix C has to be such that all the eigenvalues of CG have 
negative real parts. In order to understand what properties of G determine the choice 
of a particular stabilising transformation C, we use the following parametrisation for 
the general two-dimensional orthogonal matrices: 

-, \ ri ^ ( ^ '"'^^ ^ 

(^■t J Gg a. — 

^ ' ' \ — ssma cosa 

where s = ±1 and — tt < a < tt. When a — — 7r/2, 0, 7r/2, or tt, we obtain the set of 
matrices Csd- For example, 7r/2 = (_°io) ^^'^ C'-i.tt = (gi*!)- 

If we write G = g^, = 1, 2), then the eigenvalues of Gg^aG are given by the 
following equations: 



(2.2) cri,2 = -Acos(a - 6) ± cos^ [a - 6) - sdetG 



where det G = 511522 - 512521, A = \yj{sgii + 522)^ + (5512 - 521)^, and 

(2.3) tan0=^^i^^^, -tt < < ^ . 

-S511 - 522 

Note that the signs of the numerator and denominator are significant for defining angle 
in the specified range and should not be canceled out. It is clear from Ea. H2.2|l that 
both eigenvalues have negative real parts when 



(2.4) s = s = sgn det G, and \a - 0\ < 



2 



provided that det G 7^ 0. This result proves the validity of Conjecture ^ for n — 2. 
Moreover, it shows that there are typically two matrices in Csd that stabilise a given 
fixed point. 

Parameter 9 clearly plays an important role in the above analysis. The following 
theorems show its relationship to the eigenvalues and eigenvectors of the stability 
matrix of a fixed point. 

Theorem 2.1. Let x* be a saddle fixed point of fP{x) : i-^ whose stability 
matrix DfP(x*) has eigenvalues Ai^2 such that IA2I < 1 < |Ai| and eigenvectors defined 
by the polar angles < (/)i.2 < tt, i.e. Vi^2 = (cos 0i_2, sin (/)i_2)'''. Then the following 
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is true for angle 6 defined in Eq. I^a.yp : 
Case 1: Ai < —1 

(2.5) 

Moreover, if |Ai| S> 1, then 

(2.6) 0«(0i-02)(mod^)-f 
Case 2: Ai > 1 



(2.7) 



3ir 



2 v^l - 02 , < 01 - 02 < TT , 
f - 01 - 02 , -TT < 01 - 02 < . 



Proof. Matrix G = DfP{x*) — I, where / is the identity matrix, can be written as 
follows: 

<7ll 5l2 A _ f cos 01 cos 02 \ / Ai — 1 \ / cos 01 cos 02^ ^ 



(2 SlG ■ 

\92i 922 J \^ sin 01 sin 02 y V ^2 ~ ^ ) \ySin0i sin 02 

Case 1 : Since det G = (Ai — 1)(A2 — 1) > we set s — 1 and obtain from Eq. H2.3|l : 

(Ai - A2) COt(0i - 02) 



(2.9) tan 6* = 



2 - Ai - A2 



where, just like in Eq. (|2.3() . as well as in Eqs. (|2.10|l and (|2.13() below, the signs of 
the numerator and denominator should not be canceled out. Since 2 — Ai — A2 > 0, 
we have that cos 6* > or 



For [All > 1, Eq. (E3 yields: 

tan 6* « — cot(0i — 02) , 

and, given Eq. H2.5|) . the result in Eq. H2.6|l immediately follows. 

Case 2 : In this case detG = (Ai - 1)(A2 - 1) < 0, so, from Eq. H2.3() with s = -1: 

, . „ (A2 - Ai) cos(0i + 02)/ sin(0i - 02) 

(^•IL)) ta.nO—— , , . , , 

(A2 - Ai) sm(0i + 02)/ sin(0i - 02) 

_ - COS(01 + 02)/ Sin(0i - 02) 

- sin(0i + 02)/ sin(0i - 02) ' 

since A2 - Ai < 0. The resuh in Eq. |(22|) follows. □ 

Theorem 2.2. Let x* be a spiral fixed point of f^{x) : i-^ whose stability 
matrix DfP{x*) has eigenvalues Ai,2 = A ± iw. Then 

(2.11) ^e(-f,f) if A<1, 

e (-7r,-f) U (f,7r) if A>1. 
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Proof. The stability matrix can be decomposed as follows: 

, . nfPt *\ (coscf) ( A tj\/cos0 e'' 

(2.12) WO=(^,in^ \){s\ncp 

where 77 e R. Given that G = DfP{x*) - I, we have from Eq. 

— cj cosh 77/ sine/) 



(2.13) tan0 = 



1-A 



The result in Eq. (|2.11|) follows from the sign of the denominator. □ 

The key message of the above theorems is that the stabilising transformation ma- 
trix depends mostly on the directions of the eigenvectors and the signs of the unstable^ 
eigenvalues (or their real parts), and only marginally on the actual magnitudes of the 
eigenvalues. This means that a transformation that stabilises a given fixed point x* of 
fP will also stabilise fixed points of all periods with similar directions of eigenvectors 
and signs of the unstable eigenvalues. In the next Section, we will show how this 
observation can be used to construct stabilising transformations for efficient detection 
of periodic orbits in systems with n > 2. 

3. Extension to higher-dimensional systems. To extend the analysis of the 
preceding Section to higher-dimensional systems, we note that the matrix Cs,e, as 
defined by Eqs. ^2.1\i . (|2.3|l . and 1)2.4(1 . is closely related to the orthogonal part of 
the polar decomposition of G |1()| . Recall that any non-singular n x n matrix can be 
uniquely represented as a product 

(3.1) G = QB, 

where Q is an orthogonal matrix and _B is a symmetric positive definite matrix. The 
following theorem provides the link between Gg^g and Q for n — 2: 

Theorem 3.1. Let G gM.'^^'^ be a non-singular matrix with the polar decomposi- 
tion G — QB, where Q is an orthogonal matrix and B is a symmetric positive definite 
matrix. Then matrix Gs,e, as defined by Eqs. 1^2.'^) and i2.4\ ), is related to Q 

as follows: 

(3.2) Cs,e = -0" 



Proof. Since Gs,e is an orthogonal matrix by definition, it is sufficient to prove 
that GsfiG is symmetric negative definite. Then, by the uniqueness of the polar 
decomposition, it must be equal to —B. 

Denote by 6ij the element in the i-th row and j-th column of GgfiG. We must 
show that 612 = 621- Using Eq. 1(2. 3|l . we have that 



(3.3) 



bi2 — sgi2 cos 6 -\- g22 sin 9 
sgi2 ~ g2i 



Sgi2 + 322 

511512 



■S511 
521522 



S511 + 522 



- 522 
cos 6 



^That is, eigenvalues whose magnitude is larger than one. 



COS ( 
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and similarly 

(3.4) 621 = 521 cos — sgii sin 6* 

sgi2 - 521 



521 - sgii- 

-S.911 - 522 
<?11<?12 + 521522 



COS ( 



•S511 + 522 



cos6' , 



hence the matrix CgfiG is symmetric. Since, by definition, and s are chosen such 
that the eigenvalues of CsfiG are negative, the matrix Cs,eG is negative definite. 
Finally, by the uniqueness of the polar decomposition, 

CsfiG = —B = —Q^G, 

which completes the proof. □ 

For n > 2, we can always use the polar decomposition to construct a trans- 
formation that will stabilise a given fixed point. Indeed, if a fixed point x* of an 
n-dimensional flow has a non-singular matrix G = Dg{x*), then we can calculate the 
polar decomposition G = QB and use 

(3.5) G = -Q'^ , 

to stabilise x*. Moreover, by analogy with the two-dimensional case, we can expect 
that the same matrix G will also stabilise fixed points x with the matrix G = Dg{x), 
as long as the orthogonal part Q of the polar decomposition G ~ QB is sufficiently 
close to Q. More precisely, 

Observation 2. G will stabilise x, if all eigenvalues of the product QQ^ have 
positive real parts. 

We base this observation on the following corollary of Lyapunov's stability theo- 
rem [O] : 

Corollary 3.2. Let B G M"^" be a positive definite symmetric matrix. If 
Q G M"^" is an orthogonal matrix such that all its eigenvalues have positive real 
parts, then all the eigenvalues of the product QB have positive real parts as well. 

Proof. According to Lyapunov's theorem, a matrix A e R"^" has all eigenvalues 
with positive real parts if and only if there exists a symmetric positive definite G G 
^nxn gm^j^ i-jjj^^ ^Tq = h is positivc definite. 

Let A = QB and let's choose G in the form G = ^QB~^Q^ . Since B is positive 
definite, its inverse B~^ is also positive definite, and, since G and B~^ are related 
by a congruence transformation, according to Sylvester's inertia law |12|. G is also 
positive definite. Now, 

GQB + {QBYG = ^QB-^Q^QB + h^BQ'^QB-^Q'^ = i[Q + Q^] . 

Therefore, QB has eigenvalues with positive real parts if and only if ^[Q + Q^] is 
positive definite. The proof is completed by observing that, for orthogonal matrices, 
the eigenvalues of -I- Q^] are equal to the real parts of the eigenvalues of Q. □ 

Note that Observation |21 is a direct generalisation of conditions in Eq. (|2.4() which 
are equivalent to requiring that the eigenvalues of Gs,aCj g have positive real parts. 

In the scheme where already detected periodic orbits are used as seeds to detect 
other orbits ,4, , we can use G in Eq. (|3.5|) as a stabilising matrix for the seed x* . Based 
on the analysis in ^ this will allow us to locate a periodic orbit in the neighbourhood 
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of X* with similar invariant directions and the same signs of the unstable eigenvalues. 
Note, however, that the neighbourhood of the seed x* can also contain periodic orbits 
with the similar invariant directions but with some eigenvalues having the opposite 
sign (i.e. orbits with and without reflections). To construct transformations that would 
stabilise such periodic orbits, we can determine the eigenvalues and eigenvectors of 
the stability matrix of x* 

(3.6) DfP{x*) = VKV-^ , 

where A = diag(Ai, . . . , A„) is the diagonal matrix of eigenvalues of DfP{x*) and V is 
the matrix of eigenvectors, and then calculate the polar decomposition of the matrix 

(3.7) G = V{SA-I)V-\ 

where S — diag(±l, ±1, . . . , ±1). Note that, as follows from the analysis in [2|for 
n — 2 and numerical evidence for n > 2, changing the sign of a stable eigenvalue 
will not result in a substantially different stabilising transformation. Therefore, we 
restrict our attention to the following subset of S: 

(3.8) for . = i,...,.. 

For a seed with k real unstable eigenvalues, this results in 2*^ possible transformations. 
Note that, on the one hand, this set is much smaller than Csd, while, on the other 
hand, it allows us to target all possible types of periodic orbits that have invariant 
directions similar to those at the seed. 

4. Numerical results. In this Section we illustrate the performance of the new 
stabilising transformations on a four-dimensional kicked double rotor map |24| and a 
six-dimensional system of three coupled Henon maps [23) . Both systems are highly 
chaotic and the number of UPOs is expected to grow rapidly with increasing period. 
The goal is to locate all UPOs of increasingly larger period. Of course, the com- 
pleteness of the set of orbits for each period cannot be guaranteed, but it can be 
established with high degree of certainty by using the plausibility criteria outlined in 
the Introduction. 

In order to start the detection process, we need to have a small set of periodic 
orbits (of period p > 1) that can be used as seeds. Such orbits can be located using, 
for example, random seeds and the standard Newton-Raphson method (or the scheme 
in Eq. H1.3|l with (3 = 0). We can then use these periodic orbits as seeds to construct 
the stabilising transformations and detect more UPOs with higher efficiency. The 
process can be iterated until we find no more orbits of a given period. In our previous 
work 131 El we showed that for two-dimensional maps such as Henon and Ikeda it is 
sufficient to use period-(p — 1) orbits as seeds to locate plausibly all period-p orbits. 
For higher-dimensional systems, such as those considered in the present work, these 
seeds may not be sufficient. However, it is always possible to use more seeds by, for 
example, locating some of the period-(p + 1) orbits, which can then be used as seeds 
to complete the detection of period-p orbits. The following recipe can be used as a 
general guideline for developing a specific detection scheme for a given system: 

1. Find a set of orbit points of low period using random seeds and the iterative 
scheme in Eq. with (3 = (i.e. the Newton-Raphson scheme). 

2. To locate period-p orbits, first use period-{p — 1) orbits as seeds. For each 
seed a;o; construct 2^ stabilising transformations C using Eqs. i)^.615'. js]) . where k is 
the number of unstable eigenvalues of Df^^^ixo). 
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3. Starting from Xq and with a fixed value of f3 > use the iteration scheme in 
Eq. fl.yp to construct a sequence {xi} for each of the 2^ stabilising transformations. 
If a sequence converges to a point, check whether it is a new period-p orbit point, and 
if so, proceed to find a complete orbit by iterating the map f . 

4. Repeat steps 2 — 3 for several (3 in order to determine the optimal value of 
this parameter (see explanation below). 

5. Repeat steps 2 — 4 using newly found period-p points as seeds to search for 
period-{p + 1) orbits. 

6. Repeat steps 2 — 4 using incomplete set of period- (p + 1) orbits as seeds to 
find any missing period-p orbits. 

Although we know that the action of f3 is to increase the basin size of the stabihsed 
points, it is not known a priori what values of (3 to use for a given system and period. 
Monitoring the fraction of seeds that converge to periodic orbits, we observe that 
it grows with increasing (3 until it reaches saturation, indicating that the iterative 
scheme faithfully follows the flow S. On the other hand, larger /? translates into 
smaller integration steps and, therefore, longer iteration sequences. Thus the optimal 
value of P is just before the saturation point. As demonstrated in our previous work 0] 
and observed in the numerical examples presented in the following sections, this value 
appears to scale exponentially with the period and can be estimated based on the 
information about the detection pattern at lower periods. 

The stopping criteria in step 3, which we use in the numerical examples discussed 
below, are as follows. The search for UPOs is conducted within a rectangular region 
containing a chaotic invariant set. The sequence {xi} is terminated if (i) Xi leaves 
the region, (ii) i becomes larger than a pre-defined maximum number of iterations 
(we use i > 100 + 5/3 ), (iii) the sequence converges, such that ||5(a;i)|| < Tolg. In 
cases (i) and (ii) a new sequence is generated from a different seed and/or with a 
different stabilising matrix. In case (iii) five Newton iterations are applied to Xi to 
allow convergence to a fixed point to within the round-off error. A point x* for which 
||g(a;*)|| is the smallest is identified with a fixed point of f^. The maximum round-off 
error over the set Xp of all detected period-p orbit points 

(4.1) emax(p) = max{||5(a;*)|| : x* G Xp} 

is monitored in order to assess the accuracy of the detected orbits. 

To check if the newly detected orbit is different from those already detected, its 
distance to other orbit points is calculated: if — j/*||oo > Tolx for all previously 
detected orbit points y*, then x* is a new orbit point. Even for a large number 
of already detected UPOs, this check can be done very quickly by pre-sorting the 
detected orbit points along one of the system coordinates and performing a binary 
search for the points within Tol^ of x* . The infinity norm in the above expression is 
used for the computational efficiency of this check. 

The minimum distance between orbit points 

(4.2) d„ii„(p) =mm{\\x* -y*\\^ : x*,y* G Xp} 

is monitored and the algorithm is capable of locating all isolated UPOs of a given 
period p as long as emax(p) < Tolg < Tolx < c?min(p)- Since typically emaxip) increases 
and (imin(p) decreases with p (see Tables ETI and IT!^ . the above conditions can be 
satisfied up to some period, after which higher-precision arithmetics needs to be used 
in the evaluation of the map. For the numerical examples presented in the following 
sections we use double-precision computation with Tolg = 10^^ and Tol^ ~ 10^^. 
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Table 4.1 

Number n{p) of prime period-p UPOs, and the number N{p) of fixed points of p-times iterated 
map for the kicked double rotor map. The asterisk for p = 8 indicates that this set of orbits is not 
complete. Parameters emax(p) and dinin(p) are defined in Eqs. 14- 1( and 



p 


n{p) 


N{p) 




dmin(p) 


1 


12 


12 


1.0 ■ 10-1" 


1.3 ■ 10" 


2 


45 


102 


5.9 ■ 10-1"' 


3.4- 10-1 


3 


152 


468 


5.8 ■ 10-13 


6.2 • 10-2 


4 


522 


2190 


2.7- 10-12 


6.9 • 10-3 


5 


2200 


11012 


2.6 ■ 10-11 


1.1 ■ 10-3 


6 


9824 


59 502 


1.6 ■ 10-1° 


1.8 ■ lO-'i 


7 


46 900 


328 312 


9.7- 10-10 


9.1 • 10-5 


8* 


229 082 


1 834 566 


1.2 ■ 10-* 


5.5 • 10-5 



4.1. Kicked double rotor map. The kicked double rotor map describes the 
dynamics of a mechanical system known as the double rotor under the action of a 
periodic kick j24j . It is a four-dimensional map defined by 

o^ [ ^n+i \ ^ f Myri + Xn (mod 27r) \ 

\y7i+i ) \Lyn + csinxn+i /' 

where Xn £ are the angle coordinates and y„ £ are the angular velocities after 
each kick. Parameters L and M are constant 2x2 matrices that depend on the 
masses, lengths of rotor arms, and friction at the pivots, while c G is a constant 
vector whose magnitude is proportional to the kicking strength /q. In our numerical 
tests we have used the same parameters as in with the kicking strength /o — 8.0. 

The following example illustrates the stabilising properties of the transformations 
constructed on the basis of periodic orbits. Let us take a typical period-3 orbit point 
X* = (0.6767947,5.8315697), y* = (0.9723920,-7.9998313) as a seed for locating 
period-4 orbits. The Jacobian matrix Df^{x* ,y*) of the seed has eigenvalues A = 
diag(206.48, -13.102, -0.000373, 0.000122). Therefore, based on the scheme discussed 
in f|2lEqs. H3.6I3.8|I . we can construct four stabilising transformations C corresponding 
to {811,822) in Eq. (jSHJ) being equal to (+,+), (-,+), (+, -) and (-,-). Of the 
total of 2190 orbit points of period-4 (see Table H!T|l . the transformations Ci, C2, C3, 
and C4 stabilise #(1) = 532, #(2) = 544, #(3) = 474, and #(4) = 516 orbit points, 
respectively, and these sets of orbits are almost completely non-overlapping. That 
is, the number of orbits stabilised by both Ci and C2 is #(1 n 2) = 2. Similarly, 
#(1 n 3) = 16, #(1 n 4) = 0, #(2 n 3) = 0, #(2 n 4) = 14, and #(3 n 4) = 0. On 
the other hand, the number of period-4 orbits stabilised by at least one of the four 
transformations is #(1 U 2 U 3 U 4) = 2034. This is a typical picture for other seeds of 
period-3 as well as other periods. 

This example provides evidence for the validity of our approach to constructing 
the stabilising transformations in high-dimensional systems based on periodic orbits. 
It also shows that, in the case of the double rotor map, a single seed is sufficient for 
constructing transformations that stabilise majority of the UPOs. Of course, in order 
to locate the UPOs, we need to ensure that the seeds are in the convergence basins 
of the stabilised periodic orbits. That is why we need to use more seeds to locate 
plausibly all periodic orbits of a given period. Still, because of the enlarged basins 
of the stabilised orbits, the number of seeds is much smaller than that required with 
iterative schemes that do not use the stabilising transformations. 

Compared to the total of 384 matrices in Csd7 we use only two or four transfor- 
mations for each seed, depending on the number of unstable directions of the seed 
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orbit points. Yet, the application of the detection scheme outhned in 2] aUows us 
to locate plausibly all periodic orbits of the double rotor map up to period 7. Table 
14.11 also includes the number of detected period-8 orbits that were used as seeds to 
complete the detection of period 7. 

The confidence with which we claim to have plausibly complete sets of periodic 
orbits for each period is enhanced by the symmetry consideration. That is, since the 
double rotor map is invariant under the change of variables {x, y) ^-^ (27r — x, —y), a 
necessary condition for the completeness of the set of orbits for each period is that 
for any orbit point (a;*, y*) the set also contains an orbit point (27r — x* , —y*)- Even 
though this condition was not used in the detection scheme, we find that the detected 
sets of orbits (apart from period 8) satisfy this symmetry condition. Of course, this 
condition is not sufficient to prove the completeness of the detected sets of UPOs, but, 
combined with the exhaustive search procedure presented above, provides a strong 
indication of the completeness. 

4.2. Coupled Henon maps. Another system we use to test the efficacy of our 
approach is a six-dimensional system of three coupled Henon maps (CHM), 

(4.4) xi^^^^a-Hi^ + hx^^_^, for i = l,2,3, 

where a — 1.4 and h ~ 0.3 are the standard parameter values of the Henon map and 
the coupling is given by 

(4.5) i^„-(l-e)< + ie(a;f;^i +<-'), 

with x'~^ = x^ and = x^. We have chosen the coupling parameter e = 0.15. Our 
choice of this system is motivated by the work of Politi and Torcini \Z'6\ in which 
they locate periodic orbits in CHM for a small coupling parameter by extending the 
method of Biham and Wenzel Pp. This makes the CHM an excellent test system, 
since we can compare our results against those for the Biham- Wenzel (BW) method. 
The BW method defines the following artificial dynamics 

(4.6) xiit) ^ (-l)^("'^){x^„+i(t) - a + [S:i{t)f - bxi_,{t)]. 

with s{n,j) G {0, 1}. Given the boundary condition x^j^^ ~ x\, the equilibrium states 
of Eq. (|4.6() are the period-p orbits for the CHM. The BW method is based on the 
property that every equilibrium state of Eq. (|4.6|l can be made stable by one of the 
2'^^ possible sequences of s{n,j) and, therefore, can be located by simply integrating 
Eq. H4.6|l to convergence starting from the same initial condition x{^ — 0.0. It is 
also found that, for the vast majority of orbits, each orbit is stabilised by a unique 
sequence of s(n, j). 

In order to reduce the computational effort Politi and Torcini suggest reducing 
the search to only those sequences s{n,j) which are allowed in the uncoupled system, 
i.e. with e = 0. This reduction is possible because the introduction of coupling has 
the effect of pruning some of the orbits found in the uncoupled Henon map without 
creating any new orbits. 

We have implemented the BW method with both the full search and the reduced 
search (BW-r) up to as high a period as is computationally feasible (see Table ^2J. 
In the case of the full search we detect UPOs up to period 8 and in the case of the 
reduced search up to period 12. The seed x{^ — 0.0 was used for all periods except 
for period 4, where it was found that with this seed both BW and BW-r located only 
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Table 4.2 

The number of prime UPOs for the system of three coupled Henon maps ( CHM) detected by 
three different methods: BW - full Biham- Wenzel, B W-r - reduced Biham- Wenzel, ST - our method 
based on stabilising transformations, Max - maximum number of detected UPOs obtained from all 
three methods and the system symmetry. See text for details. 



p 


BW 


BW-r 


ST 


Max 


emax(p) 




n(p) 


1 


8 


8 


8 


8 


1.3 ■ 10"^* 


9.9 


10-i 


2 


28 


28 


28 


28 


4.6 ■ 10-1* 


5.2 


10-1 


3 
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34 


34 


40 


40 


2.7- 10-8 


4.2 


10-2 


5 




















6 


74 


74 


72 


74 


9.5 ■ 10-^° 


8.6 


10-^ 


7 


28 


28 


28 


28 


1.0 ■ 10-8 


5.6 


10-^ 


8 


271 


271 


285 


286 


1.1 • 10-'' 


5.5 


10-3 


9 




63 


64 


66 


9.9 • lO-'^ 


2.6 


10-1 


10 




565 


563 


568 


1.3 ■ 10-8 


4.1 


10-* 


11 




272 


277 


278 


7.1 ■ 10-^ 


5.4 


lo-* 


12 




1972 


1999 


1999 


2.5 • 10-'' 


4.3 


lo-* 


13* 






1079 




8.6 • 10-8 


4.0 


10-* 


14* 






6599 




2.3 • 10-'^ 


3.5 


10-* 


15* 






5899 




7.0 • 10-6 


1.5 


10-* 



28 orbits. We found a maximum of 34 orbits using the seed = 0.5. It is possible 
that more orbits can be found with different seeds for other periods as well, but we 
have not investigated this. The example of period 4 illustrates that, unlike for a single 
Henon map, the Biham- Wenzel method fails to detect all orbits from a single seed. 

Even though our approach (labeled "ST" in Table is general and does not 
rely on the special structure of the Henon map, its efficiency far surpasses the full 
BW method and is comparable to the reduced BW method. Except for periods 6 and 
10, the ST method locates the same or larger number of orbits.^ 

Unlike the double rotor map, the CHM possesses very few periodic orbits for small 
p, particularly for odd values of p. Therefore, we found that the direct application of 
the detection strategy outlined at the beginning of f^l would not allow us to complete 
the detection of even period orbits. Therefore, for even periods p we also used p + 2 
as seeds and, in case of period 12, a few remaining orbits were located with seeds of 
period 15. We did not attempt to locate a maximum possible number of UPOs for 
p > 12. The numbers of such orbits (labeled with asterisks) are listed in Table for 
completeness. 

As with the double rotor map, we used the symmetry of the CHM to test the 
completeness of the detected sets of orbits. It is clear from the definition of the CHM 
that all its UPOs are related by the permutation symmetry (i.e., six permutations 
of indices j). The column labeled "Max" in Table lists the maximum number of 
UPOs that we were able to find using all three methods and applying the permutation 
symmetry to find any UPOs that might have been missed. As can be seen in Table^21 
only a few orbits remained undetected by the ST method. 

Concluding this Section, we would like to point out that the high efficiency of 
the proposed method is primarily due to the fact that each stabilising transformation 
constructed based on the stability properties of the seed orbit substantially increases 
the basins of convergence of orbits stabilised by this transformation. This is apparent 

^The precise reason for the failure of the ST method to detect all period 6 and 10 orbits needs 
further investigation. We believe that the orbits that were not detected have uncharacteristically 
small convergence basins with any of the applied stabilising transformations. 
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in a typical increase of the fraction of converged seeds with the increasing value of 
parameter f3 in Eq. For example, when detecting period-10 orbits of CHM using 

period- 12 orbits as seeds, the fraction of seeds that converge to periodic orbits grows 
from 25-30% for small (3 (essentially the Newton-Raphson method) to about 70% for 
the optimal value of /3. 

5. Discussion and Conclusions. We have presented a new scheme for con- 
structing stabilising transformations which can be used to locate periodic orbits in 
chaotic maps with the iterative scheme given by Eq. (|1.3|l . The scheme is based on 
the understanding of the relationship between the stabilising transformations and the 
properties of eigenvalues and eigenvectors of the stability matrices of the periodic or- 
bits. Of particular significance is the observation that only the unstable eigenvalues 
are important for determining the stabilising transformations. Therefore, unlike the 
original set of transformations proposed by Schmclchcr and Diakonos, which grows 
with the system dimensionality as 2"7i!, our set has the size of at most 2'^, where k 
is the maximum number of unstable eigenvalues (i.e. the maximum dimension of the 
unstable manifold). It is also apparent that, while the SD set contains a large frac- 
tion of transformations that do not stabilise any UPOs of a given system, all of our 
transformations stabilise a significant subset of UPOs. The dependence of the num- 
ber of transformations on the dimensionality of the unstable manifold rather than 
on the system dimensionality is especially important in cases when we study low- 
dimensional chaotic dynamics embedded in a high-dimensional phase space. This is 
often the case in systems obtained from time-space discretisation of nonlinear partial 
differential equations (e.g. the Kuramoto-Sivashinsky equation). Application of the 
stabilising transformations approach to such high-dimensional chaotic systems will be 
the subject of our future work. 

The new transformations were tested on two systems: a kicked double rotor 
map and three symmetrically coupled Henon maps. We aimed to achieve a plausibly 
complete detection of periodic orbits of low periods up to as high a period as was 
computationally feasible. In both cases our algorithm was able to detect large numbers 
of UPOs with high degree of certainty that the sets of UPOs for each period were 
complete. We have used the symmetry of the systems in order to test the completeness 
of the detected sets. On the other hand, when the aim is to detect as many UPOs 
as possible without verifying the completeness, the symmetry of the system could be 
used to increase the efficiency of the detection of UPOs: once an orbit is detected, 
additional orbits can be located by applying the symmetry transformations. 

One apparent drawback of the new scheme is that a small set of UPOs needs 
to be available for the construction of the stabilising transformation at the start of 
the detection process. With the systems studied so far, we had no problem detecting 
UPOs of low period using the standard Newton-Raphson method by setting {3 = 
in Eq. 1)1. 3|l . However, in systems where it is hard to detect even a single periodic 
orbit, it would be useful to be able to determine stabilising transformations without 
the knowledge of UPOs. Since the stabilising transformations depend mostly on the 
properties of the unstable subspace, and since the decomposition into stable and un- 
stable subspaces can be defined at any, not just periodic, point on the chaotic set, it 
should be possible to estimate such properties and construct the stabilising transfor- 
mations without the knowledge of the UPOs. The decomposition could be done, for 
example, in a process similar to that used in the subspace iteration algorithm [17j . 
and a set of stabilising transformations, for example Csd, could then be applied only 
within the unstable subspace. The feasibility of such a construction will be the topic 



14 



J. J. CROFTS AND R. L. DAVIDCHACK 



for future investigation. 
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